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Abstract 

The  main  goal  of  this  joint  project  between  the  Naval  Postgraduate  School  and  the  University 
of  California,  Santa  Cruz,  was  to  analyze  the  mathematical  proprieties  of  pseudospectral  (PS) 
optimal  control  theory.  The  secondary  goal  of  our  project  was  to  exploit  the  new  properties  to 
develop  or  enhance  computationally  efficient  control  algorithms.  These  goals  were  motivated 
by  a  need  to  generate  real-time  solutions  to  constrained  nonlinear  control  problems  that  plays 
a  central  role  in  low-cost  operations  of  high-performance  space  systems.  In  this  project,  we 
developed  a  general  framework  for  PS  optimal  control  theory  based  on  discretization  over  an 
arbitrary  (i.e.  non-Gaussian)  grid.  Our  analysis  on  dual  consistency  revealed  some  hidden 
properties  of  PS  methods;  and  led  to  a  new  PS  scheme  utilizing  primal-only  conditions.  These 
new  developments  enhance  the  algorithmic  efficiency  that  is  crucial  for  real-time  applications. 
During  this  research  period,  we  also  got  an  opportunity  to  apply  our  ideas  to  the  time-optimal 
reorientation  of  a  NASA  spacecraft  in  orbit;  and  performed  a  historic  command  and  control  of 
the  Transition  Region  and  Coronal  Explorer  (TRACE)  space  telescope. 


Introduction 


High  performance  of  space  systems  demands  full  considerations  of  the  nonlinearity  and  strin¬ 
gent  constraints.  Optimal  control  plays  a  central  role  in  designing  efficient  feedback  control 
algorithms  for  such  complicated  high-performance  space  systems.  As  a  means  to  address  the 
challenges  of  future  space  systems,  we  address  the  problem  of  generating  a  real-time  solution  to 
the  following  generic  constrained  nonlinear  optimal  control  problem: 


Minimize  J[.x(-),  «(•)]  = 


Subject  to 


x(t)  = 
e(x(-l),x(l))  = 
h(x(t. ),  u(t.))  < 


E(x(- l),x(l))  +  /  F(x(t),u(t))  dt 

f(x(t),u(t )) 

0 

0 


where  F  :  RNx  x  RNu  ^  M,  E  :  RNx  x  RNx  ^  M,  /  :  RNx  x  RNu  ^  RNx,  e  :  RNx  x  RNx  i->  RNe, 
h  :  RNx  X  RNw  i— >■  RNh,  and  N  with  the  appropriate  subscript  is  some  given  natural  number. 

Over  the  last  decade,  pseudospectral  (PS)  methods  have  emerged  as  one  of  the  most  efficient 
numerical  methods  for  solving  constrained  nonlinear  optimal  control  problems,  as  evident  from 
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its  current  widespread  use,  include  experimental  demonstrations  and  flight  operations.  Based 
on  different  orthogonal  polynomials  and  Gaussian  grids,  various  PS  optimal  control  methods 
can  be  developed.  In  PS  methods,  the  grid  is  typically  determined  by  the  roots  or  the  ex¬ 
trema  of  orthogonal  functions.  Since  different  orthogonal  polynomials  have  different  properties 
with  respect  to  interpolation,  integration,  and  discretization  of  the  optimal  control  problem, 
the  approximations  of  the  states  and  costates  based  on  these  various  grids  could  have  different 
outcomes.  For  example,  Chebyshev  PS  method  offers  rapid  node  computation  and  derivative 
computation  via  a  fast  Chebyshev  differentiation  scheme  that  is  similar  to  an  FFT  computa¬ 
tion.  However,  the  non-unity  weight  function  associated  with  Chebyshev  polynomials  creates 
difficulties  in  constructing  a  proper  dual  space.  As  a  mean  to  analyze  the  effects  of  different 
orthogonal  polynomials  and  grids  on  solving  optimal  control  problems,  we  focus  on  developing 
PS  methods  over  arbitrary  grids  for  Problem  B.  Such  research  can  provides  a  unified  framework 
for  the  analysis  of  different  PS  methods.  It  also  provides  a  way  to  compare  performances  among 
different  PS  methods  and  suggest  guidelines  for  choosing  the  proper  grids  and  discretization 
approaches  for  constrained  optimal  control  problems. 


Major  Accomplishments  of  The  Project 


PS  Method  over  Arbitrary  Grids 

Because  Problem  B  has  boundary  conditions,  we  choose  to  =  —1  and  tjv  =  1  to  be  the  boundary 
points.  We  further  choose  an  arbitrary  grid  between  the  two  end  points,  and  denote  it  as 
—  l<ti<t-2<---<  tjv-i  <  1-  In  Ref. [12] ,  we  developed  a  PS  discretization  of  Problem  B  that 
generates  a  sequence  of  finite  dimensional  optimization  problems  given  by, 


limmize 


JN(x,u )  =  Y^k=oF(xk’uk)wk  +  E(x°,XN) 


Subject  to  ]C=0  DkjX3  =  f(xk, 


u 


e(x°,xN)  =  0 
h(xk,  uk)  <  0 

k  =  0, 1, . . . ,  N 


In  Problem  BN ,  D  is  a  square  differentiation  matrix  and  Wk  are  quadrature  weights  defined  by 

Dij  =  wk  =  [  4>k(t)dt 
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where  4>j(t)  are  the  Lagrange  polynomials  given  as 

9N(t ) 


N 


-  tj) 


9N(t )  =  IJ(t-ij). 


3= 0 


Note  that  both  Wk  and  D  are  constants  that  depend  on  the  grid  points  only.  Therefore,  the 
performance  of  the  PS  methods  is  determined  by  the  choice  of  the  grid.  This  general  formulation 
helps  us  to  analyze  the  effects  of  grid  selection.  For  example,  the  feasibility  and  the  consistency 
of  PS  methods  over  arbitrary  grids  can  be  guaranteed  under  certain  conditions.  One  of  them 
is  that  the  weights,  Wk,  k  =  0, 1, . . . ,  N,  need  to  be  positive  for  all  k.  When  a  uniform  grid  is 
chosen,  when  N  >  11,  at  least  one  of  the  weights  is  negative.  Therefore,  uniform  grid  is  not  an 
appropriate  choice,  as  convergence  is  not  guaranteed.  This  point  is  demonstrated  by  way  of  a 
linear  optimal  control  problem  in  Ref.  [12], 
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Dual  Consistency 

For  Problem  BN  to  converge  to  Problem  B,  the  discretization  scheme  must  commute  with  the 
dualization  process.  That  is  the  KKT  conditions  of  Problem  BN  must  preserve  the  structure 
of  the  continuous  necessary  conditions.  For  dual  consistency  analysis,  consider  a  distilled  scalar 
unconstrained  optimal  control  problem  with  no  cost  function.  The  costate,  A (t),  must  satisfy 
the  adjoint  equation  A  =  ~  ^yjy^A,  whose  PS  discretization  over  arbitrary  grids  is 

DX  =  —A(x,u)-X  (1) 

where  A  =  [Ao,---  ,Xn]t  and  A(.x,u)  is  a  diagonal  matrix  with  the  ith  diagonal  entries  be 
On  the  other  hand,  the  discretized  state  equation  Dx  =  f(x,u )  generates  the  KKT 

conditions 

Dt X  =  A(x,  u)  ■  X  (2) 

where  A  are  KKT  multipliers.  Clearly,  A  is  not  the  discretization  of  the  continuous  costate, 
since  DT  is  not  a  differentiation  operator.  However,  it  can  be  shown  that  differentiation  matrix 
defined  on  arbitrary  grid  can  always  be  factorized  into  D  =  M~1S  with  AT  and  S  be  square 
matrices  defined  by 

N  N 

Mi,j  =  sitj  = 

k= 0  fc=0 

where  ri  are  Legendre-Gauss-Lobatto  (LGL)  nodes  and  uA  are  LGL  weights.  Matrix  M  is 
always  positive  definite,  and  satisfies  Sylvester-like  equation 

MD  +  DtM  =  A  (3) 

with  all  entries  of  A  be  zero  except  for  Aoo  =  —  1  and  A^vw  =  1-  With  the  help  of  (3),  the  KKT 
system  (2)  is  transferred  into 

D  •  (A/-1  A)  =  -  [AT1  A(x,  u)M\  ■  (Af_1A)  +  Af_1A  •  (Af_1A)  (4) 

Since  the  discrete  differentiation  is  applied  to  Af_1A,  the  scaled  KKT  multiplier,  Af_1A,  serves 
as  the  discretization  of  the  continuous  costate.  This  is  the  first  hidden  property  of  PS  optimal 
control  methods:  for  PS  methods  over  arbitrary  grid  distribution,  costate  is  formed  by  the  mul¬ 
tiplication  of  adjoint  weight  matrix  M  and  the  KKT  multipliers.  For  standard  LGL  PS  method, 
since  A I  is  a  diagonal  matrix  formed  by  LGL  weights,  we  recover  the  standard  approach  for 
costate  computation,  i.e. ,  scale  the  KKT  multipliers  by  LGL  weights. 

From  (4),  it  is  also  clear  that  to  maintain  the  structure  of  the  vector  field  in  the  adjoint 
equation,  we  need  Af-1A(x,  u)M  =  A (x,u).  For  nonlinear  problems,  this  condition  holds  if  and 
only  if  A I  is  diagonal.  It  is  another  hidden  property  of  PS  methods  been  discovered.  From  the 
definition  of  M  and  the  property  of  Lagrange  interpolating  polynomial,  A 1  in  general  cannot 
be  diagonal  for  arbitrary  grid.  It  is  diagonal  only  if  the  grid,  tk,  coincide  with  LGL  points,  rk, 
i.e.,  LGL  grid  is  adopted.  It  explains  the  superior  performance  of  LGL  PS  method.  In  fact, 
the  diagonal  property  of  M  is  deteriorated  as  the  grid  distribution  is  further  away  from  the  set 
of  LGL  points.  Shown  in  Fig.l  are  the  absolute  value  of  off-diagonal  entries  of  matrix  M  with 
LGL  grid,  CGL  grid  and  uniform  grid.  It  is  clear  that  for  CGL  grid,  off-diagonal  entries  are 
small  implying  that  M  is  approximately  diagonal.  For  uniform  grid,  A I  is  far  away  from  being 
diagonal.  Since  all  Gaussian  quadrature  grids  share  similar  distribution,  such  Gaussian  grids 
maintain  the  diagonal  property  of  A 1  approximately.  Our  results  also  explain  the  connection 
between  the  Clenshaw-Curtis  weights  and  the  covector  mapping  theorem  for  the  Chebyshev  PS 
method  developed  in  Ref.  [2],  That  is,  for  a  CGL  grid,  M  is  not  diagonal  but  is  well-approximated 
by  a  diagonal  matrix  composed  of  the  Clenshaw-Curtis  weights. 
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Primal-only  Closure  Conditions 

The  last  term  in  (4)  relates  to  the  transversality  conditions.  For  LGL  PS  method,  utilizing  the 
diagonal  property  of  M,  a  set  of  primal-dual  closure  conditions  are  introduced  to  imposed  right 
transversality  condition.  This  approach  was  generalized  to  Chebyshev  PS  methods  with  M  be 
approximated  by  a  diagonal  matrix  composed  of  the  Clenshaw-Curtis  weights.  For  arbitrary  grid 
PS  methods,  to  handle  the  last  term  in  (4)  and  the  transversality  conditions,  we  modify  Problem 
Bn  as  follows:  introduce  new  primal  decision  variables  x°  and  xN ,  which  represent  almost  the 
same  quantities  as  x°  and  xN .  In  the  discretization  of  the  event  cost,  E(x(—  l),x(l)),  and  event 
conditions,  e(— 1, 1)  =  0,  we  use  x°  and  xN  instead  of  x°  and  xN .  Therefore,  the  discrete  cost 
and  the  event  conditions  are 

N 

JN  =  J2F(xk,uk)wk  +  E(x°,xN) 

k= 0 

0  =  e(x°,xN) 


The  discretization  of  the  dynamical  equations  is  modified  as 

Dx  =  f(x,u)  +  (x°  —  x°)ro  +  (xN  —  xN)rN 


where  f(x,u)  =  [f(x°,  u°),  ■  ■  ■  ,f(xN,uN)]T.  ro  and  rjy  are  the  first  and  the  last  column  of 
M^1.  It  can  be  shown  that  ro  and  r n  are  explicitly  given  by 


ro,k 

rN,k 


(— l)^1  . 

“ — 2 - LN(tk)(t  -  1) 

—LN(tk)(t  +  1),  k  =  0, 1,  •  •  ■  ,  N 


where  L]\r(t)  is  the  Nth  order  Legendre  polynomial.  Such  modification  generates  a  modified 
Problem  BN  that  enables  dualization  to  commute  with  discretization.  This  new  scheme  facili¬ 
tates  enhancing  the  algorithmic  efficiency  of  the  underlying  spectral  algorithm  by  generating  a 
full-rank  KKT  system.  The  proposed  primal-only  closure  condition  is  easy  to  implement  and 
applicable  to  different  PS  methods  including  Legendre  and  Chebyshev  PS  methods  as  special 
cases. 


Flight  Results 

During  this  research  period,  we  also  got  an  opportunity  to  apply  our  ideas  to  the  time-optimal 
reorientation  of  a  NASA  spacecraft  in  orbit.  On  August  2010,  we  led  a  team  of  engineers  and 
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operators  to  perform  a  historic  command  and  control  of  the  Transition  Region  and  Coronal 
Explorer  (TRACE)  space  telescope.  Prior  to  our  flight  test,  no  spacecraft  had  ever  performed 
a  minimum-time  attitude  maneuver.  This  is,  in  part,  because  of  the  well-known  difficulties  in 
solving  the  minimum-time  spacecraft  attitude  control  problem.  Our  techniques  allow  an  operator 
to  perform  verification  and  validation  of  the  command  signals  without  a  need  to  understand  PS 
optimal  control  theory.  Because  TRACE  is  actuated  by  reaction  wheels,  the  resulting  dynamics 
are  not  only  non-Eulerian  but  also,  the  control  space  is  non-convex  and  state-dependent.  The 
well-known  difficulties  in  solving  such  problems  were  easily  circumvented  by  our  PS  theory.  The 
mathematical  implications  as  well  as  the  engineering  results  acquired  from  the  telemetry  are 
being  analyzed  now  with  the  goal  of  publications  in  the  coming  years. 

Other  Advancements 

During  this  research  period,  we  also  extended  our  previous  work  ou  feasibility,  consistency  and 
convergence  analysis.  When  the  Legendre  PS  method  is  used  for  feedback  linearizable  systems, 
we  showed  in  Ref.  [3]  that  the  rate  of  convergence  is  at  least  ^2/73-1  ,  where  m  is  defined  by  the 
smoothness  of  the  optimal  trajectory;  i.e.  an  appropriate  Sobolev  space.  If  the  optimal  control 
is  C00,  then  the  convergence  rate  can  be  made  faster  than  any  given  polynomial  rate.  Relative 
to  our  earlier  results,  we  also  weakened  the  assumptions  in  the  theorems  we  proved.  In  this 
research  project,  we  also  analyzed  the  performance  of  a  PS  algorithm  for  generating  real-time 
solution  of  infinite-horizon  optimal  control  problems.  We  proved  the  closed-loop  stability  under 
box  constraints  on  states  and  controls.  Such  results  are  important  for  real-time  space  control 
applications. 


Ongoing  Work  and  Future  Research  Plans 

This  is  the  final  year  of  a  three- year  project.  The  new  developed  PS  computational  scheme  paves 
the  way  to  construct  more  efficient  algorithms  for  solving  optimal  control  problems,  for  example, 
multiscale  PS  methods  for  dynamical  systems  with  different  timescales.  Many  space  applications 
involve  both  slow  and  fast  dynamics.  The  nature  of  such  systems  with  mixed  time-constants 
suggests  a  possible  efficient  discretization  using  multi-resolution  for  different  subsystems.  The 
new  developed  arbitrary  grid  PS  method  provides  flexibility  in  designing  different  grids  for 
different  subsystems.  We  will  explore  this  idea,  together  with  a  covector  mapping  theorem  for 
costate  computation,  to  develop  multiscale  PS  methods. 
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